
#include <stdlib.h>
#include <stdio.h>

#include "globals.h"
#include "geometry.h"
#include "polygon.h"
#include "distance.h"
#include "sgidrawi.h"

#ifdef Graphics
#include "gl/gl.h"
#include "gl/device.h"
#include "forms.h"
#endif

static int      printInit = 0;

/*	Allocate memory space for a new polygon.  The sequence of vertices
	are copied to the polygon data structure and the 'vertices' 
	parameter can be used for something else. */
Pgon 
*newPolygon(int numPoints, Pt *vertices)
{
	Pgon	*poly;
	Pt		*pt;
	int		i;
	Coordinate	extrL, extrR, extrH, extrB;
	Coordinate	myX, myY;
	
	poly = (Pgon *) malloc (sizeof(Pgon));
	pt = (Pt *) malloc (sizeof(Pt) * numPoints);

	if ((poly != NULL) && (pt != NULL)) {
		poly->points = pt;
		poly->numPoints = numPoints;

		/*	Set up default bounding box values. */

		poly->leftmost = 0;
		poly->rightmost = 0;
		poly->upper = 0;
		poly->lower = 0;
		poly->voronoiCount = 0;

		extrL = XVal(vertices[0]);
		extrR = extrL;
		extrH = YVal(vertices[0]);
		extrB = extrH;

		/*	Copy all the points to the polygon structure, watching for
			the maximum and minimum value entries. */

		for (i = 0; i < numPoints; i++) {
			PointAssign(poly->points[i], vertices[i]);
			myX = XVal(vertices[i]);
			myY = YVal(vertices[i]);

			if (myX > extrR) {
				poly->rightmost = i; extrR = myX;
			} 
			else if ((myX < extrL) || 
			((myX == extrL) && (myY > YVal(vertices[poly->leftmost])))) {
				poly->leftmost = i; extrL = myX;
			}

			if (myY > extrH) {
				poly->upper = i; extrH = myY;
			} 
			else if (myY < extrB) {
				poly->lower = i; extrB = myY;
			}
		}
	}
	else if (poly != NULL)
		free (poly);

	return (poly);
}

/*	Find the middle vertex of the polygon given two indices of the
	polygon.  The procedure assumes, for now, that the initial indices
	provided are valid for the polygon. */
int	
middlePt (Pgon *poly, int lower, int upper)
{
	return ((upper - lower) % poly->numPoints);
}

/*	Release the memory associated with a polygon */
void 
killPolygon(Pgon **polygon)
{
	if (*polygon != NULL) {
		(*polygon)->numPoints = 0;	/* In case there is a handle lying about */
		free ((*polygon)->points);
		free (*polygon);
		*polygon = NULL;
	}
}

/*  These routines are just for printing the polygons */

/*	Allocate display resources to draw polygons.  This includes
	creating a window for the polygons. */
void 
initPolyPrint(void)
{
    if (printInit == 0) {
#ifdef Graphics
        fl_init();
        fl_set_graphics_mode(FALSE, FALSE);
        create_free_forms();
		clearScreen();
		swapbuffers();
#endif
        printInit = 1;
    }
}

/*	Release all display resources reserved for drawing polygons. */
void 
endPolyPrint(void)
{
    if (printInit == 1) {
#ifdef Graphics
		clearScreen();
		swapbuffers();
        winclose(doodle);
#endif
        printInit = 0;
    }
}

/*	Procedure to draw a polygon on the terminal screen.  The routine
	initPolyPrint must have been invoked before this point */
void 
printPoly(int colour, Pgon *polygon)
{
#ifdef Graphics
	int	i;

	if (printInit != 0) {
		if (polygon->numPoints > 1) {
			for (i = 1; i < polygon->numPoints; i++) {
				drawline (colour, (float)XVal(polygon->points[i-1]),
					(float)YVal(polygon->points[i-1]),
					(float)XVal(polygon->points[i]),
					(float)YVal(polygon->points[i]));
			}
			drawline (colour, 
				(float)XVal(polygon->points[polygon->numPoints - 1]),
				(float)YVal(polygon->points[polygon->numPoints - 1]),
				(float)XVal(polygon->points[0]),
				(float)YVal(polygon->points[0]));
		}
		else {
			drawpoint (colour, (int) XVal(polygon->points[0]),
					(int) YVal(polygon->points[0]));
		}
	}
#endif
}



/*	Possible answers:
		0:  point right of both lines.
		1:	point left of ln1, right of ln2.
		2:	point right of ln1, left of ln2.
		3:	point left of both lines.
*/

int 
pointConfig(Line *ln1, Line *ln2, Pt *testPt)
{
	int	ans;

	ans = 0;

	if (LineSide(ln1, testPt) == LeftSide)
		ans += 1;
	if (LineSide(ln2, testPt) == LeftSide)
		ans += 2;

	return (ans);
}

#define modPt(A,B) ((A)->points[(B)%(A)->numPoints])

void 
nearestPoints (Pgon *p, Pgon *q, Pt *ppt, Pt *qpt, int *pidx, int *qidx)
{
	Line	pccwNorm, pcwNorm, qccwNorm, qcwNorm, tmpLine;
	int		phigh, plow, pmid, qhigh, qlow, qmid, result;
	Pt		ignoreme;

	outerCommonTangent (p, q, &ignoreme, &ignoreme, &plow, &qhigh);
	outerCommonTangent (q, p, &ignoreme, &ignoreme, &qlow, &phigh);
	if (plow >= phigh)
		phigh += p->numPoints;
	if (qlow >= qhigh)
		qhigh += q->numPoints;

	while ( (phigh - plow > 1) ||
			(qhigh - qlow > 1) ) {
		pmid = (plow + phigh) / 2;
		qmid = (qlow + qhigh) / 2;

		/*	Get the normals at pmid. */

		LineOf(&(modPt(p, pmid-1 + p->numPoints)), &(modPt(p, pmid)), &tmpLine);
		perpRightTurnLine(&(modPt(p, pmid)), &tmpLine, &pccwNorm);

		LineOf(&(modPt(p, pmid)), &(modPt(p, pmid+1)),&tmpLine);
		perpRightTurnLine(&(modPt(p, pmid)), &tmpLine, &pcwNorm);

		/*	Get the normals at qmid. */

		LineOf(&(modPt(q, qmid-1 + q->numPoints)), &(modPt(q, qmid)), &tmpLine);
		perpRightTurnLine(&(modPt(q, qmid)), &tmpLine, &qccwNorm);

		LineOf(&(modPt(q, qmid)), &(modPt(q, qmid+1)),&tmpLine);
		perpRightTurnLine(&(modPt(q, qmid)), &tmpLine, &qcwNorm);

		/*	Narrow down the binary search interval. */

		result = pointConfig(&pccwNorm, &pcwNorm, &(modPt(q, qmid)))*4 + 
			pointConfig(&qccwNorm, &qcwNorm, &(modPt(p, pmid)));
		switch (result) {
			case 0:
			case 8:
			case 4:
				plow = pmid;
				break;
			case 11:
			case 15:
			case 5:
			case 6:
			case 7:
				phigh = pmid;
				break;
			case 2:
			case 3:
			case 1:
			case 9:
				qlow = qmid;
				break;
			case 12:
			case 14:
			case 13:
				qhigh = qmid;
				break;
			case 10:
				plow = pmid; phigh = pmid;
				qlow = qmid; qhigh = qmid;
				break;
			default:
				printf ("Reached an unexpected case in finding nearest "
						"points of polygons %d.\n", result);
				plow = pmid; phigh = pmid;
				qlow = qmid; qhigh = qmid;
				break;
		}
	}

	/*	We're down to two line segments to handle in the worst case. */

	phigh = phigh % p->numPoints;
	plow = plow % p->numPoints;
	qhigh = qhigh % q->numPoints;
	qlow = qlow % q->numPoints;

	if (phigh == plow) {
		ptNearestSegment (&(p->points[plow]), &(q->points[qlow]),
			&(q->points[qhigh]), qpt);
		PointAssignPtr(ppt, &(p->points[plow]));
		*pidx = plow;
		*qidx = qlow;
	} 
	else if (qhigh == qlow) {
		ptNearestSegment (&(q->points[qlow]), &(p->points[plow]),
			&(p->points[phigh]), ppt);
		PointAssignPtr(qpt, &(q->points[qlow]));
		*pidx = plow;
		*qidx = qlow;
	}
	else {
		Real	dists[4], maxdist;
		Pt		actual[4], near[4], start[2], end[2];	
		int		idx[4], i, maxidx, ignore;

		/*	Really brute force for now. */
	
		PointAssign(actual[0], p->points[plow]);
		PointAssign(actual[1], p->points[phigh]);
		PointAssign(actual[2], q->points[qlow]);
		PointAssign(actual[3], q->points[qhigh]);

		PointAssign(start[0], q->points[qlow]);
		PointAssign(start[1], p->points[plow]);
		PointAssign(end[0], q->points[qhigh]);
		PointAssign(end[1], p->points[phigh]);

		for (i = 0; i < 4; i++) {
			ptNearestSegment (&(actual[i]), &(start[i/2]),
				&(end[i/2]), &(near[i]));
			squarePointDistance(&(actual[i]), &(near[i]), &(dists[i]));
		}

		maxidx = -1;
		for (i = 3; i >= 0; i--) {
			if ((maxidx == -1) || (dists[i] < maxdist)) {
				maxidx = i;
				maxdist = dists[i];
			}
		}

		if (maxidx == -1) {
			printf ("Trouble -- nobody close enough to take the blame\n");
		}
		else {
			if (maxidx < 2) {
				PointAssignPtr(ppt, &(actual[maxidx]));
				PointAssignPtr(qpt, &(near[maxidx]));
			}
			else {
				PointAssignPtr(qpt, &(actual[maxidx]));
				PointAssignPtr(ppt, &(near[maxidx]));
			}
		}

		/*	Sluff these off for now. */

		*pidx = plow;
		*qidx = qlow;
	}

}


